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New Developments in Relativistic Viscous Hydrodynamics 

Paul Romatschke^ 

Institute for Nuclear Theory, University of Washington, Box 351550, Seattle, WA, 98195 

Starting with a brief introduction into the basics of relativistic fluid dynamics, 
I discuss our current knowledge of a relativistic theory of fluid dynamics in the 
presence of (mostly shear) viscosity. Derivations based on the generalized second law 
of thermodynamics, kinetic theory, and a complete second-order gradient expansion 
are reviewed. The resulting fluid dynamic equations are shown to be consistent for all 
these derivations, when properly accounting for the respective region of applicability, 
and can be applied to both weakly and strongly coupled systems. In its modern 
formulation, relativistic viscous hydrodynamics can directly be solved numerically. 
This has been useful for the problem of ultrarelativistic heavy-ion collisions, and I 
will review the setup and results of a hydrodynamic description of experimental data 
for this case. 
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I. INTRODUCTION 

A. Non-relativistic fluid dynamics 

Fluid dynamics is one of the oldest and most successful theories in modern physics. In its 

non-rclativistic form, it is intuitively understandable due to our everyday experience with 
hydrodynamics, or the dynamics of water^. The degrees of freedom for an ideal, neutral, 
uncharged, one-component fluid are the fluid velocity v{t,x), the pressure p{t,x), and the 
fluid mass density p{t,x), which are linked by the fluid dynamic equations [1],[2]§2, 

dtv+(v-d)v — —-dp, (1) 

dtp + pd ■ V + V ■ dp = . (2) 

These equations are referred to as "Euler equation" (1) and "Continuity equation" (2), 
respectively, and typically have to be supplemented by an equation of state p = p{p) to close 
the system. For non-ideal fluids, where dissipation can occur, the Euler equation generalizes 
to the "Navier-Stokes equation" [3, 4],[2]§15, 

dv'^ J. dtf 1 dp 1 (9n*^' . . 

dt dx^ p dx^ p dx^ ' 

^\dx'^^ dxi 3 dxi) ^ dx^ ' ^ ' 

where Latin indices denote the three space directions, e.g. i = 1,2,3. The viscous stress 
tensor 11^* contains the coefficients for shear viscosity, r], and bulk viscosity, C, which are 
independent of velocity. The non-relativistic Navier-Stokes equation is well tested and found 
to be reliable in many applications, so any successful theory of relativistic viscous hydrody- 
namics should reduce to it in the appropriate limit. 



B. Relativistic ideal fluid dynamics 

For a relativistic system, the mass density p{t, x) is not a good degree of freedom because 
it does not account for kinetic energy that may become sizable for motions close to the 
speed of light. Instead, it is useful to replace it by the total energy density e(t,x), which 
reduces to p in the non-relativistic limit. Similarly, v(t,x) is not a good degree of freedom 
because it does not transform appropriately under Lorentz transforms. Therefore, it should 
be replaced by the Lorentz 4-vector for the velocity, 

where Greek indices denote Minkowski 4-space, e.g. p = 0,1,2,3 with metric g^,j = 
diag(-|-, — , — , — ) (the same symbol for the metric will also be used for curved spacetimes). 



^ In some fields it has been the tradition to use the term hydrodynamics synonymous with fluid dynamics 
of other substances, and I will adopt this somewhat sloppy terminology. 
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The proper time increment dT is given by the hne element, 

{dTf = g^.^dx^dx'' = {dtf - {dxf , 

2" 



dx 



{dtf 1 - {vf 



where here and in the following, natural units h — c — — 1 will be used. This implies 
that 

' (6) 



dT dt yjl-v"^ 



which reduces to = {l,v) in the non-relativistic limit. In particular, one has = (1,0) 
if the fluid is locally at rest ("local rest frame"). Note that the 4- vector only contains 
three independent components since it obeys the relation 



2 u V 



(7) 



so one does not need additional equations when trading v for the fluid 4- velocity u^. 

To obtain the relativistic fluid dynamic equations, it is sufficient to derive the energy- 
momentum tensor T^*^ for a relativistic fluid, as will be shown below. The energy-momentum 
tensor of an ideal relativistic fluid (denoted as T^q^) has to be built out of the hydrodynamic 
degrees of freedom, namely two Lorentz scalars (e,p) and one vector u^, as well as the 
metric tensor g^^j. Since T^^ should be symmetric and transform as a tensor under Lorentz 
transformations, the most general form allowed by symmetry is therefore 



rj) = e {cog"'' + ci^'^w'^) +p{c2g'''' + c^u'^u^) 



(8) 



In the local restframe, one requires the T^Q^j component to represent the energy density e 
of the fluid. Similarly, in the local rest frame, the momentum density should be vanishing 
^(0) ^ 0' space-like components should be proportional to the pressure, T^q^ = p5^^ 

[2] §133. Imposing these conditions onto the general form (8) leads to the equations 



(co + ci)e (c2 + C3)p = e. 



-coe - C2P = p. 



(9) 



which imply cq = 0, ci = 1, C2 = -1, C3 = 1, or Tj^ = eui'vy - p (g^"" - u^u''). For later 
convenience, it is useful to introduce the tensor 
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HI/ 



(10) 



It has the properties A^^u^ = A^^'^u^ = and A^'^A" = A^" and serves as a projection 
operator on the space orthogonal to the fluid velocity u'^. In this notation, the energy- 
momentum tensor of an ideal relativistic fluid becomes 



^'^^euV-pA^r 
If there are no external sources, the energy-momentum tensor is conserved, 

f) — n 
^M-' (0) - • 



(11) 



(12) 



5 



It is useful project these equations in the direction parallel (u^dnT^) and perpendicular 
(A"9^T(qP to the fluid velocity. For the first projection, one finds 

= {e + p)dX + n''d^e = Q , (13) 

where the identity u^dijU^ — \dfj_{ui,u'^) — \dijl — was used. For the other projection one 
finds 

= (e + p) u^'d^.u'^ - A^'^di.p = . (14) 

Introducing the shorthand notations 

D = ui'd^, V" = A^"9^ (15) 

for the projection of derivatives parallel and perpendicular to u^, equations (13), (14) can be 
written as 

De + {e + p)^^u^' = (16) 
(e+p)^^" - V> = . (17) 

These arc the fundamental equations for a relativistic ideal fluid. Their meaning becomes 
transparent in the non-relativistic limit: for small velocities 1^1 <C 1 one finds 

D^u^'di.^dt + v-d + Oilvl^), ^ A'^'di.c^d' + 0{\v\), (18) 

so D and V* essentially reduce to time and space derivatives, respectively. Imposing further 
a non-relativistic equation of state where p ^ e, and that energy density is dominated 
by mass density e ~ p, Eq. (16) becomes the continuity equation (2), and Eq. (17) the 
non-relativistic Eulcr equation (1). 

One thus recognizes the fluid dynamic equations (both relativistic and non-relativistic) 
to be identical to the conservation equations for the fluid's energy-momentum tensor. 



II. RELATIVISTIC VISCOUS HYDRODYNAMICS 

A. The relativistic Navier-Stokes equation 

In the ideal fluid picture, all dissipative (viscous) effects arc by definition neglected. If 
one is interested in a fluid description that includes for instance the effects of viscosity, one 
has to go beyond the ideal fluid limit, and in particular the fluid's energy momentum tensor 
will no longer have the form Eq. (11). Instead, one writes 

where is the familiar ideal fluid part given by Eq. (11) and 11'^'' is the viscous stress 
tensor that includes the contributions to T'^'^ stemming from dissipation. Considering for 
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simplicity a system without conserved charges (or at zero chemical potential), all momentum 
density is due to the flow of energy density 

u^T^"" ^eu"" — > u^n^'' = 0. (20) 

While here this is the only possibility, for a more general system with conserved charges 
one can view this as a choice of frame for the definition of the fiuid 4-velocity, sometimes 
referred to as Landau-Lifshitz frame. This can be easily understood by recognizing that in 
a system with a conserved charge there will be an associated charge current that can be 
used alternatively to define the fiuid velocity, e.g. the Eckart frame u^n^ — n. These choices 
refiect the freedom of defining the local rest frame either as the frame where the energy 
density (Landau-Lifshitz) or the charge density (Eckart) is at rest. Since the physics must 
be the same in either of these frames, one can show that charge diffusion in one frame is 
related to heat fiow in the other frame, as done e.g. in the appendix of [5]. For other recent 
discussions of relativistic viscous hydrodynamics in the presence of conserved charges, see 
e.g. [6, 7]. 

Similar to the case of ideal fiuid dynamics studied in section IB, the fundamental equa- 
tions of viscous fiuid dynamics are found by taking the appropriate projections of the con- 
servation equations of the energy momentum tensor, 

u.d^T^''' = De + {e + p)dX + w.<9^n^" = , 

A'^d^T"'' = (e + P)^«" - V> + A'^d^U'"' = . (21) 

The first equation can be further simplified by rewriting Ui.d^Il'^'^ — 8^ [u^Jl^^) — W^^d(^ptUi,), 
and using the identity 

= u^D + (22) 

as well as the choice of frame, m^II^'^ = 0. Here and in the following the (. . .) denote 
symmetrization, e.g. 

Hence, the fundamental equations for relativistic viscous fiuid dynamics are 



De + {e + p)d^u^-WV(^u,) = 0, 

(e + p)Du'' - V> + Ay^W = . (23) 



At this point, however, the viscous stress tensor has not been specified. Indeed, much of 
the remainder of this work will deal with deriving expressions for li^^ , which together with 
(23) will give different theories of viscous hydrodynamics. 

An elegant way of obtaining H'^'^ builds upon the second law of thermodynamics, which 
states that entropy must always increase locally. The entropy density s is connected to energy 
density, pressure and temperature T by the basic equilibrium thermodynamic relations for 
a system without conserved charges (or zero chemical potential), 

e + p = Ts, Tds = de . (24) 

The second law of thermodynamics can be recast in the covariant form 



daS^ > 



(25) 
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using the entropy 4-current s'^ which in equihbrium is given by 

s'^ = su^" . (26) 
The thermodynamic relations (24) allow to rewrite the second law (25) as 

d^s^^ = Ds + sd^u^ = ^De + ^-^d^u^ = ^n^'^V(^«,) > , (27) 

where (23) was used to rewrite De. It is customary to split 11'*'' into a part tt'"' that is 
traceless, tt^ = 0, and a remainder with non-vanishing trace, 

j^nv ^ ^ . (28) 

Similarly one introduces a new notation for the traceless part of V(jj,u,^), 

V<^ii,> = 2V(^M,) - ^A^,V«m" , (29) 

so that the the second law becomes 

d^s^" = ^7r''''V<^ii.> + ^nV„ii" > . (30) 

One recognizes that this inequality is guaranteed to be fulfilled if 

n"'' ^ rjV^^u''^ , U^CVaU"', rj > , C>0, (31) 

because then d^s^ is a positive sum of squares. 

In the non-relativistic hmit, the viscous stress tensor becomes that of the Navier-Stokes 
equations (4), which leads one to equate rj,C. with the shear and bulk viscosity coefficient, 
respectively. Also, for this reason we refer to the system of equations (23), (28), (31) as the 
relativistic Navier-Stokes equation. While beautifully simple, it turns out that the relativistic 
Navier-Stokes equation - unhke its non-relativistic counterpart - exhibits pathologies for all 
but the simplest flow profiles, as will be shown below. 

B. Acausality problem of the relativistic Navier-Stokes equation 

Let us consider small perturbations of the energy density and fiuid velocity in a system 
that is initially in equilibrium and at rest, 

e^eo + Se{t,x), ^ (1,0) + Su''{t,x), (32) 

where for simplicity the perturbation was assumed to be dependent on one space coordinate 
only. The relativistic Navier-Stokes equation then specifies the space-time evolution of the 
perturbations. For the particular direction a — y, Eq. (23) gives 

(e + p)Duy - V^p + Aldf.n^" = (eo + po)dtSuy + d^U'^ + 0{6^) , 
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This implies a diffusion-type evolution equation for the perturbation 5u''^{t,x): 

dt5vy "^dlSu^ = 0{5^) . (33) 

Co + Po 

To investigate the individual modes of this diffusion process, one can insert a mixed Laplace- 
Fourier wave ansatz 

into Eq. (33). This gives the "dispersion- relation" of the diffusion equation, 

u = -^^k' , (34) 

which one can use to estimate the speed of diffusion of a mode with wavenumber k, 

v^(k) = ^ = 2^i^k . (35) 
^ ^ dk eo+po ^ ' 

One finds that vt is linearly dependent on the wavenumber, which implies that as k becomes 
larger and larger, the diffusion speed will grow without bound. In particular, at some 
sufficiently large value of k, VT{k) will exceed the speed of light, which violates causality^. 
Therefore the relativistic Navier-Stokes equation does not constitute a causal theory. 

The obvious conclusion to draw from this argument is that the relativistic Navier-Stokes 
equation exhibits unphysical behavior for the short wavelength (/c 3> 1) modes and hence 
can only be valid in the description of the long wavelength modes. This is not a principal 
problem, as one can regard hydrodynamics simply as an effective theory of matter in the 
long wavelength, k ^ limit. However, having a finite range of validity in k typically 
is a practical problem when dealing with more complicated flow profiles that do not lend 
themselves to analytic solutions and have to be solved numerically. In this case, it turns out 
that the high k modes are associated with instabilities [8] that make it necessary to regulate 
the theory by other means. A simple argument to understand the practical problem can 
be given as follows: modes that travel faster than the speed of light in one Lorentz frame 
correspond to modes traveling backwards in time in a different frame. Hydrodynamics is 
an initial value problem which requires a well defined set of initial conditions. However, if 
there arc modes present in the equations that travel backwards in time, the initial conditions 
cannot be given freely [9], and as a consequence one cannot solve the relativistic Navier- 
Stokes equation numerically. 

One possible way to regulate the theory is provided by considering the "Maxwell-Cattaneo 
law" [10, 11] 

r^dtU-y + U-y^-r)odJuy (36) 

instead of the Navier-Stokes equation. Here is a new transport coefficient sometimes 
referred to as relaxation time. The effect of this modification becomes apparent when recal- 
culating the dispersion relation for the perturbation 6u'^ using Eq. (36). One finds 

^0 k'^ 



^ One should caution that the diffusion speed exceeding the speed of hght is a hint - but no proof - of 
causality violation. The proof is given in the appendix. 
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which coincides with the dispersion relation of the diffusion equation Eq. (34) in the hydro- 
dynamic {uj, k ^ 0) hmit. More interestingly, however, is that for large frequency uj ^ 1 
Eq. (37) does not describe diffusive behavior, but instead propagating waves with a propa- 
gation speed that is finite in the limit of A; ^ 1, 

v^- ^ hm = J-^^ , (38) 

unless — > 0. Interestingly, for all known fluids the limiting value ■\J^^~^^^^ has been 
found to be smaller than one, so that the Maxwell-Cattaneo law seems to be an extension 
of the Navier- Stokes equation that preserve causality^. 

For heat flow, the corresponding Maxwell-Cattaneo law implies a dispersion relation 
equivalent to Eq. (37), and there the propagating waves can be associated with the phe- 
nomenon of second sound [12] §4, [13], observed experimentally in solid helium [14]. It is not 
known to me whether propagating high frequency shear waves, as suggested by Eq. (37), 
have been found in experiments. 

While the Maxwell-Cattaneo law seems to be a successful phenomenological extension 
of the Navier-Stokes equation, it is unsatisfactory that Eq. (36) does not follow from a 
flrst-principles framework, but is rather introduced "by hand". It will turn out, however, 
that the Maxwell-Cattaneo law - while not derivable - does seem to correctly capture some 
important aspects of relativistic viscous hydrodynamic theory, for instance that terms of 
higher order in k (higher order gradients) are needed to restore causality. 



C. Miiller-Israel- Stewart theory, entropy-wise 

In section II A, the Navier-Stokes equation was derived from the second law of thermo- 
dynamics d^s^ > by using the form of the entropy current in equilibrium, = su^. 
However, it is not guaranteed that the entropy current equals its equilibrium expression for 
a dissipative fluid that can be out of equilibrium. Speciflcally, it was suggested [15, 16] that 
out of equilibrium the entropy current can have contributions from the viscous stress tensor, 
which is sometimes referred to as "extended irreversible thermodynamics" [12, 17]. Assum- 
ing that the entropy current has to be algebraic in the hydrodynamic degrees of freedom 
and that deviations from equilibrium are not too large so that high order corrections can be 
neglected, the entropy current has to be of the form [7, 15, 16] 

= SU>^ - ^U^n' - ^U^TT^^TT'^'' + 0{U') , (39) 

where (3o,P2 are coefficients that quantify the effect of these second-order modifications of 
the entropy current. Using again Eq. (23) to rewrite d^s^ as in section II A one finds 



+ " (v^u^ - hi TD - P^Dn - ^/3on9x) > . (40) 



^ See the appendix for a proof of causality. 
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The inequality is guaranteed to be fulfilled if 



n = c fv«M^ 



T 



(41) 



with 77, C the usual bulk and shear viscosity coefficients. Note that Eq. (41) coincides with 
the Navier-Stokes equation in the limit of Po, ^2 0. For non-vanishing ^0,^2, Eq- (41) 
contains time derivatives of tt^^, fl, which are similar (but not identical) to the Maxwell- 
Cattaneo law Eq. (36) if one identifies P2 = (and similarly, (3o = ^). The set of equations 
(23), (41) (and some variations thereof) are commonly referred to as "Miiller-Israel-Stewart" 
theory and will be discussed more in section III. 

Similar to section II B, one can study the causality properties of the Miiller-Israel-Stewart 
theory by considering small perturbations around equilibrium, Eq. (32). Keeping only per- 
turbations to first order, Eq. (23) and Eq. (41) become 

dtSe + (eo + Po)dJu^ = 0, (cq + Po)dtSu^ + d.p + ^^OT^" = 0, 



(eo + Po)dtSuy + df,6Wy = 0, 5U^"' = ^tt^^^ + g^"'5U 



(42) 



The equation of state e = e{p) relates the pressure and energy density gradients, dxP = ^d^e, 
and the condition ■u^II^'^ = implies 611^" = (9(5^). Using a Fourier ansatz 



in Eq. (42) gives the system of equations 



ik 



iuj Seu,,k - ik{eo+Po) Sul ,, = 0, 
'4 ikrjo ik(o 



3 1 -|- iujTT^ 
iuj{eQ +pq) 6ul,^ - ik 



+ 



1 -|- icuTu ^ 
1 -|- iujT^ 



(43) 
(44) 

(45) 



Eq. (45) corresponds to result from the Maxwell-Cattaneo law for the transverse velocity 
perturbation Su^, discussed in section II B. The other two equations correspond to density 
perturbations and longitudinal fiuid velocity displacements, commonly known as sound. The 
sound dispersion relation is given by 



.k^ dp ,2/^4 

lou — I — + K I -- 

uj de 



and in the hydrodynamic limit {u!,k <^ 1) becomes 

'2 ^0 1 Co \ 



+ 



Co 



= 0. 



(46) 



±kcs + ik 



k^ 
'2^ 



^3eo + Po 2eo + Po, 

2 



Vo 



+ 



Co 



3eo+po 2eo+po 



2c 



Vo 



Co 



3 eo + po 2 eo + Po 



Tn 



+ Oik"") . (47) 
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The quantity 

(48) 

can be recognized to be the speed of sound when calculating the group velocity limji.-+o 
For large wavenumbers and frequencies, Eq. (46) gives a limiting sound mode group velocity 

k^oo dk V 3T^(eo+po) Tn(eo+po) 
which together with the result for the transverse mode Eq. (38) suggests that the Miiller- 
Israel-Stewart theory - derived via an extended second law of thermodynamics - constitutes 
a relativistic theory of viscous hydrodynamics that obeys causality if the relaxation times 
Tt^, Tji are not too small. Note that the requirement < 1 from Eq. (49) is more restrictive 
than Eq. (38) concerning the allowed values of c^, rj, (, r^, tu- 

However, many questions remain unanswerable within this formalism, e.g. how to obtain 
the value of Ttt, rn, or whether the assumption that the entropy current should be algebraic in 
the hydrodynamic degrees of freedom is valid (Refs. [18-20] seem to indicate the contrary). 
Therefore, it is necessary to have a different derivation of viscous hydrodynamics. 



III. FLUID DYNAMICS FROM KINETIC THEORY 

A. A very short introduction to kinetic theory 

Kinetic theory treats the evolution of the one-particle distribution function f{p,t,x), 
which can be associated with the number of on-shell particles per unit phase space, 

dN 

If coUisions between particles can be neglected, the evolution of / follows from Liouville's 
theorem, 

dT dT dt dT' dx ^ ' 

Multiplying (51) by the mass m of a particle and recognizing m-j^ = m'y{v) = — 
mvj{v) = p as the particle's energy and momentum, respectively, Eq. (51) becomes 

P^dJ = Q, (52) 

where has to fulfill the on-shell condition p'^p^ = . 

Eventually, collisions between particles cannot be neglected, and hence Eq. (51) is no 
longer valid. Taking the effect of collisions into account changes the evolution equation 
[21]§3 to 

p''d,f^-C[fl (53) 

where C[f] is the collision term that is a functional of / and the precise form of which 
depends on the particle interactions. Eq. (53) is known as the "Boltzmann-cquation" [22]. 
For a system in global equilibrium / is stationary, /(p, x) — f{o)ip) so that the Boltzmann 
equation gives 

P^5^/(o)=0 = -C[/(o)], 
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which imphes that the colhsion term vanishes in equihbrium. Note that this means that 
Eq. (52) holds for two very different regimes, namely firstly when one can ignore collisions 
(and the system is typically far from equilibrium) and secondly when the collisions are 
strong enough to keep the system in equilibrium. The first case is typically applicable if 
the timescale of the description is short enough so that the effect of particle collisions can 
be neglected. Ultimately, however, particle collisions will become important and drive the 
system towards equilibrium. It is this second case, or more generally the long time (small 
frequency, long wavelength) limit that corresponds to hydrodynamics (see also the discussion 
in section II B). 

Given the interpretation of / in Eq. (50), the particle number density should be propor- 
tional to / d^pf, or the sum of / over all momenta with weight unity. Summing instead with 
a weight of particle energy / d^p p^f, one expects a result proportional to the product of 
number density and energy, or energy density, which is a part of the energy-momentum ten- 
sor. More rigorously, one can define the relation between the particle distribution function 
and the energy-momentum tensor [23] as 

where the l.h.s. again can be understood as a sum over momenta, with the ^-function 
imposing the condition of only counting on-shell particles and the step-function to restrict 
the sum to positive energy states. 



B. Ideal fluid dynamics from kinetic theory 

In the following, I will limit myself to considering the ultrarelativistic limit where all 
particle masses can be neglected, m — 0. From Eq. (54), this leads to = 0, or vanishing 
conformal anomaly. Interpreting (54) as the fluid's energy-momentum tensor, this amounts 
to setting the bulk viscosity coefficient to zero, C = (cf. Eq. (19,31) and the discussion in 
section IV E). 

Introducing for convenience the shorthand notation 

/c?X^^5(p%)2^(/), (55) 
and taking the first moment of the Boltzmann equation, one finds 

/ dxp'p^djip'^, x^) = -j dxp'Cif] = d,J dxp'p'-fip, x) = d.T'^^ . (56) 

For particle interactions that conserve energy and momentum, the integral over the collision 
term vanishes, / dxp^C[f] = 0. If T'^" can be interpreted as a fluid's energy-momentum 
tensor, then this implies that the first moment of the Boltzmann equation corresponds to 
the fundamental equations of fluid dynamics (23), since these follow from d^T^^ = 0. 

The interpretation of the kinetic theory energy-momentum tensor in the fluid picture is 
most transparent in equilibrium, where f{p,t,x) = /(o)(p)- Similar to the discussion in the 
introduction, f(o){p) is not an optimal description for a relativistic system, since it is not 
manifestly invariant under Lorentz transformations. It is better to trade /(o) with a more 
convenient function, 

f{0){p) /eq (^^) ' 
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where u'^ is a four vector that reduces to — > (1)0) the restframe of the heat bath with 
temperature T. Eq. (54) can then be written as 

T(^oI = / ^XP>Veq if^) = a^ouV + a2i A'^^ (57) 

where in hindsight it is more convenient to choose -u^-u*^, A'^'^ as a tensor basis then u^^u", g^^ . 
The coefficients a2o, 021 are functions of the temperature only and are determined by con- 
tracting (57) with u'^u^ and A'*'', respectively, 

a2o = / dx{p^u,ff^ (^) , ^2i^\jdx {fp, - ip^u.r) /eq (^) . (58) 

Identifying with the fluid four velocity, Eq. (57) corresponds to the ideal fluid energy- 
momentum tensor Eq. (11) with e = 020, P = —021, and the equation of state e = 3 p 
(or speed of sound squared — |) following from on-shell condition in the massless limit. 

To calculate 020,021, one has to specify the equilibrium distribution function /eq. A 
concrete example where the evaluation is straightforward is for a single species of particles 
that obey Boltzmann statistics, so that /eq = exp — ■ In this case, 020 is 

easily calculated by choosing the convenient frame i*'^ = (1, 0), so that 

and 021 = —§020 — For a single species of particles obeying Bose- Einstein statistics, 

feq{x) — [e^ — 1]~^, the result would be 020 = ^^^o~- The relation between 020 and e can be 
used to re-express the temperature in terms of the energy density. 



C. Out of equilibrium 

From Eq. (57) it is evident that when the argument of the distribution function / depends 
only on scalars and one Lorentz vector u'^, the form of the energy-momentum tensor for 
kinetic theory is the same as for ideal fluid dynamics. If the system is locally in equilibrium, 
/cq is completely characterized by a vector-valued function that specifies the local rest frame 
of the heat bath, u^{x), and the local temperature (or energy density). Therefore, a system 
that is in perfect local equilibrium is described by ideal fluid dynamics. Departures from 
equilibrium result in departures from the ideal fluid dynamics picture, and hence can only 
be captured with dissipative (viscous) fluid dynamics. One can derive the correspondence 
between kinetic theory out of equilibrium and viscous hydrodynamics by considering small 
departures from equilibrium where 

/(P^ = /eq (^) [1 + Sfip^, x^] , (59) 

and 5f{p'^, x'^) -C 1. Using Eq. (59) in the definition of the energy momentum tensor Eq. (54) 
and demanding that it should correspond to Eq. (19) from viscous hydrodynamics, one finds 

j.^. = TJ^ + j dxp^p'^U^Sf = Ti;^ + 1:^^ , 

ni^^ = / dxp'p'^f^df , (60) 
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where again the contribution IT proportional to bulk viscosity was dropped because of the 
ultrarelativistic limit in Eq. (55). The out-of-equilibrium correction to the distribution 
function 6f may depend on scalars, the heat bath vector u^, the metric Qj^^, and gradients 
thereof. To make progress, it is convenient to make the dependence of Sf on the momentum 
exphcit, e.g. by using a truncated expansion in a Taylor-hke series [24] 

Sf{p^, x'') = c + p'^c^ + py'c^p + 0{p^) , (61) 

or using a different basis [23]. Using this expression in Eq. (60) and integrating over mo- 
menta, one can proceed to obtain the coefficients c, Cq, Cap in a (somewhat tedious) cal- 
culation [24]. A more direct way (that gives the same result) is to assume - similar to 

section TIC - that Sf must be an algebraic function of the hydrodynamic degrees of free- 
dom, e,p,u^ , g^'^ ,71^'^ . Then the requirement that 5f vanishes in equilibrium implies that 
c = 0, Cq = 0, and Cq,^ = C2Tiap with C2 a function of the thermodynamic variables e,p. The 
relation Eq. (60) then leads to 

tt'^'^ = 7^«/3C2/^'^"^ (62) 
where corresponds to the n — A case of the integral definition 



j dxp^'p^' ...p^-h^. (63) 



Note that for the special case of two indices I^^ = T^^, and a decomposition into Lorentz 
tensors similar to Eq. (57) can be done for each of the integrals (63). In particular, one finds 

If^'^^P = aioU^u''u''u^+a^i (u^'u^A'"^ + perm.) +042 (A'^'^A"'^ + A^^A'^^ + A/^^A^^) , (64) 

where "perm." denotes all non-trivial index permutations. Contracting the indices in 
Eq. (62) using the properties of the shear part of the viscous stress tensor, u^n^'^ — 0,7r^ — 0, 
one finds C2 = and, with Eq. (59), the distribution function for small departures from 
equilibrium takes the form 



/(p^a;-) = /e.(^) 



1 + -^p-p/^ 

2a42 



(65) 



The coefficients 040,041,042 can be calculated the same way as 020, O21 in section IIIB once 
/eq has been specified. For the special case of a Boltzmann gas where feq{x) — e~^, a 
straightforward calculation gives the relation 

042 = (e + P)T'^ , 

which holds also when allowing for nonzero particle masses. 

The equation (65) establishes the relation of the particle distribution function out of 
(but close to) equilibrium and viscous hydrodynamics. Still missing is for a relation of 
the Boltzmann equation (53) and viscous hydrodynamics is an expression for the collision 
term. Depending on the nature of the particle interactions, C[f] will have a particular, 
and sometimes complicated, form that can be simplified by assuming small departures from 
equilibrium, cf. Eq. (59). 

If one identifies the magnitude of 6f with the size of gradients of hydrodynamic degrees 
of freedom, a shortcut to obtain C[f] to lowest order in a gradient expansion is to insert (59) 
into the Boltzmann equation: 

C[f] = -p'd, [/eq (1 + 5f)] = -p^d.f^ + 0{S') . (66) 
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This approach is similar to the Chapman-Enskog approach to fluid dynamics [25]. 

For the special case of particles obeying Boltzmann statistics, /eq = e~^, the calculation of 
C[f] from Eq. (66) to first order in gradients is simple and will be given here as an illustrative 
example. Using the fundamental equations of viscous fluid dynamics (23), one can rewrite 

P'd.f.^ if^) = -P V/eq (V, + U,D) ^ 

/eq (V^«. + «[^V,i InT + ^«^«.Va«") + 0{5^), (67) 



T 

where here and in the following [. . .] denote antisymmetrization, e.g. 

Since the structure p'^p'^ in Eq. (67) is symmetric in the indices and vanishes when contracted 
with Qi^i, because of the on-shell condition, Eq. (67) reduces to 

P^a^/eq = -^V<^W.>/eq + 0{S') , (68) 

where V<^Mi/> was deflned in Eq. (29). To flrst order in gradients, the Navier-Stokes equation 
(31) is valid and as a consequence one flnds 

C[f] = |^^.^/eq + 0{S') , (69) 

which establishes the relation between the collision term and viscous hydrodynamics to flrst 
order in gradients. 



D. Miiller-Israel- Stewart theory, kinetic theory-wise 

In a theory with conserved charges where / dxC = 0, the integral over momenta (or zeroth 
moment) of the Boltzmann equation gives 

J dxp^d^f ^d^j dxp^f = d^n^' = 0, (70) 

or conservation of charge current n^. The first moment of the Boltzmann equation (shown 
in Eq. (56)) gives the conservation of the energy- momentum tensor, since / dxp°'C = 0. 
However, the integral / dxp°'p^C does not trivially vanish, unless the system is in equilibrium. 
Therefore, the second moment of the Boltzmann equation 

/ dxpVp^d^f = - / dxpy^Cif] , (71) 

will carry some information about the non-equilibrium (or viscous) dynamics of the system 
[26]. Considering again small departures from equilibrium, Eq. (65) implies 

J dxpy^dj = (/"^'^ + ^^"''^"') , (72) 
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where the integrals /^i/^a-. /in ^ere defined in Eq. (63). Similar to Eq. (64) one can do a 
tensor decomposition of J"^'^, jais^iyS^ with coefficients 030,031 and a50jO5ijO52j respectively. 
To extract the relevant terms from Eq. (72), it is useful to use a tensor projector on the part 
that is symmetric and traceless, 



(73) 



with properties u"P^^ = u^^P^^ = 0, A"^P^^ = and tt^^P^^ = 271"". Using this projector 
on Eq. (72), one finds after some algebra 



P^ad^r^'^ = PZa^i [/^A"^ + 2V("«'^)1 = 2a3iV<^M"> , 



2(342 



^52 
042 



37r("'^ 



(74) 



To calculate the r.h.s. of Eq. (71) one would have to specify the precise form of the 
collision integral. If one is only interested in the form of the equation, not the coefficients 
of the individual terms, it is convenient to again assume Boltzmann statistics, feq{x) = e~^, 
because then the form of the collision term is given by Eq. (69) and one finds 



2a427r' 



2Trj Trj Tr] 

The coefficients 031, 042, 052 are readily evaluated for a massless Boltzmann gas. 



(75) 



031 



42^5 42^6 24r^ 

^42 = TT-i 052 



and after collecting the terms from Eq. (75) and Eq. (75) one finds for the second moment 
of the Boltzmann equation (71) the result 



042 



u 



(76) 



It is useful to rewrite the expression P^^n'^^V tf,u°' in this equation by introducing the fluid 
vorticity 

^afS = V[aM/3] , (77) 

which is antisymmetric, flai3 = —^/Sa- After some algebra one finds the relation 



= -27r'^('^Ql + 



27] 6 



(78) 



where (31) was used to rewrite V^ipu^y to first order in gradients. For the massless Boltz- 
mann gas, one furthermore has DlnT = Dlne^^^ — — IV^m" + O^d"^), so that Eq. (76) 
becomes 
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= riV<V> + 0{5^) , (79) 
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where the expression was labeled to make the connection to the Maxwell-Cattaneo 

"42 

law Eq. (36) exphcit. Eq. (79) constitutes a different variant of the Miiller-Israel- Stewart 
theory, and the connection between this equation and Eq. (41), which was derived earlier in 
section II C via the second law of thermodynamics, will be discussed below. 

Since T,r multiplies all the terms in Eq. (79) which are at least of second order in gradients, 
it is a generalization of the concept of hydrodynamic transport coefficients (such as shear 
viscosity 77), and is accordingly referred to as a second order transport coefficient. For a 
Boltzmann gas, the known values of 042, 052 imply the relation 

which together with = | give the definite values v^^^ = and vf^^^ = ^ for the 
propagation speeds Eq. (38,49). This indicates that the theory by Miiller, Israel and Stewart 
does indeed preserve causality since signal propagation is subluminal. 

For a massive Boltzmann gas, one can recalculate the coefficients 052, 042 to show that 
the more general relation 

Tj e + p 

V ^ Q _L Tds_ 

s dT 

holds. Also, for Bosc-Einstcin statistics, the proportionality factor | in Eq. (80) is only 
changed by a few percent [27]. Thus it seems that the property of causality of the viscous 
fluid dynamic equations (23), (79) is fairly robust whenever kinetic theory is applicable. 



E. Discussion of Miiller-Israel- Stewart theory 

Eq. (79) contains the Navicr-Stokcs equation (31) in the limit of small departures from 
equilibrium where second order gradients (all the terms multiplied by in Eq. (79)) can be 
neglected. However, the form of the terms to second order in gradients is such that Eq. (79) 
reproduces the phenomenologically attractive feature of the Maxwell-Cattaneo law, namely 
flnite signal propagation speed. In addition, the kinetic theory derivation of Eq. (79) also 
gives a definite relation between shear viscosity and relaxation time, the first and second 
order transport coefficients, respectively, which implies not only finite, but subluminal signal 
propagation. 

However, the evolution equation for the shear stress n^'^ differs between the derivation 
from kinetic theory Eq. (79) and the second law of thermodynamics, Eq. (41), respectively. 
To make this more apparent, it is useful to rewrite Eq. (41) for the case of a Boltzmann gas 

with = ^ = 



r]V<i'vy> + 0{5'^) , (81) 



where Eq. (23) was used to rewrite DhiT = — |VaM" + C(5^). One first notes that the 
terms involving the time derivative Dtt"^ differ between Eq. (81) and Eq. (79), 

This difference is easily explained by noting that for the derivation of Eq. (81), only the 
projection of tt^j, on Eq. (81) was required to have well defined sign (40). But the differ- 
ence between Eq. (81) and Eq. (79) vanishes when contracted with 71^^;/, so these terms do 
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not actually contribute to entropy production and therefore are not naturally captured by 
the derivation in section II C. Nevertheless, one can convince oneself that these terms are 
necessary and important by contracting Eq. (79) and (81) with m^: unlike Eq. (79), the 
contraction does not vanish for (81), but instead gives Uf^Dn^" — which amounts to an 
extra (unphysical) constraint on the evolution of the shear stress tensor [28] . Therefore, the 
variant Eq. (79) of Miiller-Israel-Stewart theory derived from kinetic theory is superior to 
Eq. (81) in this respect. 

However, when inserting the kinetic theory result Eq. (79) into the conservation equation 
for the entropy current (40), one finds for the shear viscosity contribution the requirement 



2T 



>0, (82) 



where the identity n ^j,ii'^'^^VL^]p = has been used. On the one hand, there is no obvious 
reason why Eq. (82) should be fulfilled for all values of n^^, but on the other hand the second 
law of thermodynamics > should not be violated. 

One solution to the problem is that Eq. (82) may still be fulfilled if departures from 
equilibrium are small enough, so that the first term in Eq. (82) — being second order in 
gradients and manifestly positive [29] — is larger than the other term, which is third order 
in gradients, tt'^'^tt^,^ ~ 0{5'^) 0{6^). In other words, the region of apphcability of viscous 
hydrodynamics would coincide with the region of applicability of the gradient expansion 
used to derive it. 

However, most likely one also has to give up the assumption made in Eq. (39) about the 
particular form of the generalized entropy current. Indeed, a different form for s'^ allowing 
for gradients [18-20] does seem to imply d^s^ > for evolution equations of n^" that are 
more general than Eq. (81). 

While this implies that the correct theory is more complicated, Eq. (79) is a good can- 
didate for a theory of relativistic viscous hydrodynamics that fulfills the necessary minimal 
requirements, namely reduction to Navier-Stokes equation in the limit of long wavelengths, 
and causal signal propagation. The shortcoming of Eq. (79) is that the equation has un- 
known corrections to second order in gradients 0{6^), stemming from the unknown form 
of the collision term. Since the second-order gradients on the l.h.s. of Eq. (79) are needed 
to guarantee finite signal propagation, it docs not seem to be consistent to ignore terms 
of second order on the r.h.s. Rather, one would want to have a more general theory that 
includes all terms to second order in gradients. 

IV. A NEW THEORY OF RELATIVISTIC VISCOUS HYDRODYNAMICS 

A. Hydrodynamics as a gradient expansion 

All the hydrodynamic results discussed so far can be classified in terms of a gradient 
expansion of the fluid's energy momentum tensor^ T^^'^ — T^!^ + tt'*'', namely 

• Ideal Hydrodynamics: contains no gradients (zeroth order). 



* Again for simplicity only the case of shear viscosity is discussed, where 11'*'' = n^" . 
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• Navier-Stokes Equation: contains first order gradients, 

• Miiller-Israel-Stewart theory: contains second order gradients, 

As discussed in the introduction, the ideal hydrodynamic energy-momentum tensor is the 
most general structure allowed by symmetry, and therefore the zeroth order gradient expan- 
sion is complete. On the other hand, section HIE indicates that the Miillcr-Israel-Stewart 
theory potentially misses terms of second order in gradients, and hence the gradient expan- 
sion may not be complete to this order. To obtain the most general structure of viscous 
hydrodynamics to second order, one has to completely classify all possible terms in n^" to 
first and second order gradients of the hydrodynamic degrees of freedom [29]. 

To first order, since the equation of state links the pressure to the energy density, the 
only independent gradients are dfj_u"^d^je. Decomposing = + u^D^ the fundamental 
equations (23) can be used to express all time-like derivatives Du", Dt in terms of space-like 
gradients V^, so only the latter are independent. This imphes that the shear-stress tensor 
should have the structure 

TT^'^ = CM'^V"^ + CsA'^'^Va^/" + C6 U^^^V'h , (83) 

where C4, C5, Cg are functions of e only. The Landau-Lifshitz frame condition u^n'^^ = 
imphes that Cg = 0, or the absence heat flow (see section II A). Furthermore, since effects 
from bulk viscosity have been ignored, the stress tensor is traceless, which gives C5 = — |c4. 
Choosing the proportionality constant C4 = 2?], one finds ti^^ = ri'V^'^u''-^ , which shows that 
the Navier-Stokes equation corresponds to a complete gradient expansion to first order. 

To second order in gradients, the analysis proceeds similar to the one above, but there 
are more terms to consider. It turns out that for the case of only shear viscosity, there is an 
additional restriction for n^'^ besides u^tt^^^ = and tt^ = 0, namely conformal symmetry, 
that can be used to reduce the number of possible structures. 



B. Conformal viscous hydrodynamics 

A theory is said to be conformally symmetric if its action is invariant under Weyl trans- 
formations of the metric, 

g,u^g,. = e-'^^^^g,,, (84) 

where w(x) can be an arbitrary function of the spacctime coordinates, and hence g^iy is 
the metric of curved rather than flat spacetime. While on the classical level many theories 
obey this invariance, quantum correction typically spoil the symmetry, giving rise to a non- 
vanishing trace of the energy momentum tensor. One distinguishes between theories where 
in flat space quantum corrections generate 7^ — such as SU(N) gauge theories ("non- 
conformal") — and those where conformal symmetry is unbroken, such as A/" = 4 Super 
Yang- Mills ("conformal"). Note that even for "conformal" theories quantum corrections 
may couple to gravity, such that the trace of the energy-momentum tensor is non- vanishing 
in curved space ("Weyl anomaly") [30], 



(85) 
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The Weyl anomly Vr[5'^,y] in four dimensions is a function of the product of either two 
Riemann tensors R^u\p, two Ricci tensors i?^,^ or two Ricci scalars R, and hence is of fourth 
order in derivatives of g^^,, since R^^Xp, R^iv and R are all second order in derivatives [31]. 
Being interested in a gradient expansion to second order, one may therefore effectively ignore 
the presence of the Weyl anomaly. To second order in gradients, conformally invariant 
theories thus have a traceless energy-momentum tensor, which in addition transforms as 

under a Weyl rescaling in four dimensions [29] (see also the discussion in section IV E). It 
is this additional symmetry of conformal theories that helps to restrict the possible second 
order gradient terms in a theory of hydrodynamics in the presence of shear viscosity. For 
curved space, there are 8 possible contributions of second order in gradients to -k^^ that obey 

TT^^ = 0, U^Tl^^- = 0, 

DO" In e D''> In e, D<^'D''> In e, V<''u''> ( Vaii") , P^^ V<'^u^> g^s^^^u^^ 

Pap V<''u^>g^s^^\ P^^ Q'^^g^s^f^^, UjR^<^"'>^us, R^'"'^ , (87) 

but only five combinations of those transform homogeneously under Weyl rescalings, n^'^ — > 
(here and in the following Da denotes the (geometric) covariant derivative in curved 
space). The calculation is straightforward but somewhat lengthy, so I only demonstrate the 
ingredients by studying again the first order result, tt'^'^ = rjS/^^u^^ . Under conformal 
transformations, dimensionless scalars are invariant, u'^g^i,u'^ = 1 = u'^gni,u^ which implies 

under Weyl rescalings. Furthermore, the transformation of the ideal fluid's energy momen- 
tum tensor, T^^^^ = e u'^u'' -pAf"" ^ fg^ = e^^^^^Tj) then requires 

e ^ e = e^'^^^^e . (88) 

For conformal fluids, the shear viscosity coefficient is related to the energy density by 
7] oc e^/^, so that one has rj ^ fj = e'^'^^'^^^r]. Since n^'^ = rfsj^^u^^ , one then has to ver- 
ify that the first order derivative transforms homogeneously as W^^u^^ g3«;(a;)y</i^!^> 
From the expansion 

V-«'> = Wu- + Vu- - Ha-V„«" = A-0„„' + A-D^W - Ia^D^W (89) 

it becomes clear that one has to study the transformation property of the covariant derivative 
of the fluid velocity, 

where F^^ are the Christoffel symbols given by 

Ki3 = ^g"^ i9a9pi3 + dpgpa - dpg^p) ■ 

The transformation of the Christoffel is readily calculated from the transformation of the 
metric (84), 
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so that together with the transformation property of the fluid velocity one finds 

Using this result in Eq. (89) one finds that all the terms involving derivatives of the scale 
factor w{x) cancel, 

so that indeed the first order expression for tt^'^ transforms homogeneously under Weyl 

transformations. 

To second order, one repeats the above analysis for all of the eight terms in Eq. (87), 
combining them in such a way that all the derivatives of w{x) cancel. One finds the result 



TT' 



A^A^D7r"'3 + ^7r'^^(Va«") 



K 

2 



Ai 



2r)' 
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(91) 



where Ttt, k, Ai, A2, A3 are five independent second order transport coefficients, and Eq. (31) 
has been used to rewrite some expressions, disregarding correction terms of third order in 
gradients. Eq. (91) is the most general expression for tt^^ to second order in a gradient 
expansion in curved space for a conformal theory. 



C. Hydrodynamics of strongly coupled systems 

Particularly interesting examples of conformal quantum-field theories arc those that have 
known supcrgravity duals in the limit of infinitely strong coupling [32]. Since fluid dynamics 
is a gradient expansion around the equilibrium of the system, Eq. (23), (91) should be gen- 
eral enough to also capture the dynamics of these strongly coupled quantum systems in the 
hydrodynamic limit. These systems will in general not allow for a quasiparticle interpreta- 
tion, since the notion of a (quasi-)particle hinges on the presence of a well-deflned peak in 
the spectral density, which may not exist at strong coupling. Therefore, infinitely strongly 
coupled system are very different than systems described by kinetic theory (which relies on 
the presence of quasiparticles) , making their hydrodynamic limit interesting to study. 

If a known supcrgravity dual to a strongly coupled field theory is known, one can calculate 
Green's functions in these theories (for a review, see for instance [33]). A particular example 
is the Green's function for the sound mode in strongly coupled J\f = A SYM theory, with 
gravity dual on a AdS^ x S5 background, which gives rise to sound dispersion relation [29] 

k ik . 3 2 In 2,q / ■, a \ /^^\ 

uj = ±^ + — - ± — ^ k^ + Oik^) . 92 

v/3 67rT 6A/3(27rr)2 ^ ' ^ ' 

By comparing to the hydrodynamic sound dispersion relation Eq. (47), one finds the values 
for the speed of sound, shear viscosity and relaxation time for strongly coupled Af — A SYM, 

^ - ^ - ^ r - (93) 
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Calculating other quantities both in J\f — A SYM and hydrodynamics [29] and rcdcriving 
the fluid dynamic equations from the supergravity dual of A'' = 4 SYM [34], one additionally 
finds 

As a side remark, note that the dispersion relation for transverse perturbations (the shear 
mode) discussed in section II B, is ill-suited to determine the second order transport coef- 
ficients such as Tn, because information about tq enters only at fourth order in gradients 
(37), and therefore receives corrections from terms not captured by second-order conformal 
hydrodynamics [29, 35]. 

As expected, the hydrodynamic limit of strongly coupled A/" = 4 SYM reproduces the 
structure of Eq. (23), (91), which had to be true if these equations are truly universal. Fur- 
thermore, plugging the values (93) into the sound mode group velocity for large wavcmimbers 
(49), one finds ti™^^ ~ 0.92; this suggests that the hydrodynamic theory Eq. (23), (91) obeys 
causality for strongly coupled J\f = i SYM. Interestingly, this seems to be also the case 
for other known gravity duals, for instance AdSc+i, for D > 2, corresponding to strongly 
coupled conformal field theories in D spacetime dimensions. There has been an exten- 
sive amount of work on calculating the second-order transport coefficients in these theories 
[29, 34-38], which are now known analytically for all D > 2 [39] 

_ D + n[2/D - 1] . _VD - _ vn[2/D - 1] _ 

4^ ' -^^"8^' 2^ ' " - 27rT{D - 2) ' 

(95) 

with Ti[x] the harmonic number function [39, 40] 



n[x\ = / dz- = 7^; + 

Jo 1 — z 



dz 



Note that the special case D = A corresponds to the results (93) for strongly coupled 
J\f = 4 SYM, and that the ratio ^ = is universal for all of these, in line with the 
observation of Ref. [41]. Also, there seems to be some universahty for the second order 
transport coefficients: for instance, it has been found that 4Ai -|- A2 = '^V^tt for a class of 
strongly coupled field theories [42, 43]. 

Generalizing Eq. (49) to arbitrary spacetime dimension gives 



= hm 



k- 



oodk \l (D-l) r^{e + p)^ rn{e + p)- ^ ^ 



For conformal theories, C = and = and using the values (95), one finds that that 
^max jg decreasing monotonously with D > 2 from its maximum at D = 2, where y™*^^ 
would reach unity^. The values for t;™^^ and v™*^ (Eqns. (96,38), respectively) for spacetime 
dimensions D < 10 are shown in Figure 1). 

For D — A, corresponding to strongly coupled M — A SYM, also the corrections at finite 
(but large) coupling strength A to the transport coefficients have been calculated [44-46] , 



? = if (1 + f^C(3)A-3/2 + . . ^ r.T = ^ + |iC(3)A-3/2 



+ . 



^ For D = 2, the conformal field theory does not have any (first or second order) transport coefficients, but 
is completely characterized by ideal fluid dynamics [38]. 
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FIG. 1: The limiting velocities for longitudinal (96) and transverse (38) perturbations as a function 
of spacetime dimension in conformal second order hydrodynamics. 
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(l - if C(3)A-^/^ + ...), Ai = ^ (l + ^C(3)A- 



+ . . . 



(97) 



which lead to wf^^ ^ 0.92 - 0.9796A~3/2^ 



D. Hydrodynamics of weakly coupled systems and discussion 

Weakly coupled theories in general have a well defined quasiparticle structure and hence 
it is expected that the hydrodynamic properties of these theories are captured by kinetic 
theory. In particular, it is known that kinetic theory correctly reproduces the results from 
finite temperature quantum field theories, in the hard-thermal-loop (resummed one-loop) ap- 
proximation [47]. As a consequence, one would expect that the dynamics of weakly coupled 
quantum field theories in the hydrodynamic limit are well captured by the Miiller- Israel- 
Stewart theory derived via kinetic theory in section III. Comparing Eq. (91) to Eq. (79) — 
and recalhng the approximation used to derive (79) — one finds 



6 rj 



, Ai = r/r^ + T2C(l), = -2r]T^, A3 = 0, k = 



(98) 



where 0{1) refiects the unknown contribution to Ai from the collision term and k = stems 
from rederiving Eq. (79) in curved space [29]. First order transport coefficients have been 
calculated in the weak-coupling limit for high temperature gauge theories [48], in particular 
Af = 4: SYM [49] 

I ~ A2 1n(2.36A-i/2) " ^ ^ 

More recently, all second-order transport coefficients were evaluated consistently in QCD 
and scalar field theories at weak coupling [50, 51], 



5.0... 5.9 T] 



Ai 



4.1... 5.2 T]^ 

T 7 



A2 



-2r]T^, A3 = 0, K 



5 s 



(100) 
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where the range of values indicate the dependence on the couphng constant g (see Ref. [50] for 
details). Note that the results from Eq. (79) agree reasonably well with the full calculation 
Eq. (100). In particular, the value of the relaxation time is such that the limiting velocity 
^max -g sjna^gp than for strongly coupled systems. 

Comparing the results (100) to those obtained for strongly coupled theories (95), one 
finds that A3 always vanishes. This could indicate that there is an additional, unidentified 
symmetry in conformal hydrodynamics that forces this coefficient to be zero. Moreover, a 
direct calculation shows that the value of k is beyond the accuracy of kinetic theory [50-52] . 
This indicates that the kinetic theory result is not general enough to capture the dynamics of 
conformal fluids in the hydrodynamic limit for arbitrary coupling, at least when spacetime 
is curved (since k couples only to the Riemann and Ricci tensor, it does not contribute 
to Eq. (23) when spacetime is flat; however, n does enter in correlators for the energy- 
momentum tensor in flat space [29]). A possible reason for this could be the fact that 
kinetic theory itself is only a gradient expansion to first order of the underlying quantum 
field theory [47], thereby possibly missing second-order contributions. 

Furthermore, one finds that Ai^2 have the same sign for both kinetic theory and the 
strongly coupled systems studied, which could indicate that the sign of these coefficients 
does not depend on the couphng. Finally, the fact that never exceeds unity for in- 
finitely strongly coupled theories, for theories at large (but finite) coupling, and at weak 
coupling, suggests, but does not prove, that causality in a second-order conformal hydro- 
dynamics description is obeyed. At this time, there is only a proof for theories that have 
dual description in terms of Gauss-Bonnet gravity, where it has been shown that causality 
in second order hydrodynamics follows from the causality of the field theory itself [53] . As a 
consequence, one may hope that the system of equations (23), (91) constitutes a valid starting 
point to attempt a description of real (but nearly conformal) laboratory fiuids at relativistic 
speeds. This application of the viscous hydrodynamic theory to high energy nuclear physics 
will be discussed in section V. 



E. Non-conformal hydrodynamics 

Since most quantum field theories that successfully describe nature are not conformal 
theories, one may ask how deviations from conformality change the hydrodynamic equations. 
In particular, one may ask how important non-conformal terms not included in Eq. (91) are 
once conformal symmetry is slightly violated. To this end, consider the specific example of 
a SU(N) gauge theory at high temperature which has a trace anomaly [54], 

= g.^T'^^ ^ {^^g^-g-^PF^^F^^)r + W[g,^] , F;^ = d,At - d^A^ - gyu faKcA^Al , 

(101) 

where fabc are the SU(N) structure constants, ^4^ are the gauge fields and ()t denotes the 
thermal quantum field theory average. Similar to section IV B, the Weyl anomaly ^^[(7/^,^] is 
not important for what follows and will be ignored. The change of the gauge theory couphng 
gYM when changing the renormalization scale A is given by the beta-function. 



(102) 
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which for weakly coupled SU(N) gauge theories is given by [55] 

/3(^ym) = -^^ + 0(^^m)- (103) 

In fact, Eq. (101) can be derived from the gauge theory action when performing a Weyl 
transformation (84) of the partition function and noting that the renormalization scale 
changes according to A ^ e^^^^A, 

^^"^cxv/=^T-, (104) 



IJ,V 



where g is the determinant of the metric g^^, (not to be confused with gYm)- Taking another 
functional derivative of the trace anomaly [29] leads to 

= V=S [iT'Hix -y)+ sj-^) , (106) 

where the symmetry of the second derivative of the partition function with respect to the 
metric (cf. Eq. (104)) was used. On the other hand, using Eq.(lOl) one finds 



^ (v^y..(^)7^-) = ^(^^ [^FrK^ - K'Fi^) )tS{x -y)+ Oig^) , 



Sgapiy) 

= 0{g'^^) , (106) 
so that for Weyl transformations (84) this implies 



Note that an exact calculation gives = C'(fl'YM) weakly coupled SU(N) gauge theories 
[56]. Recalling that all terms in Eq. (91) transform as ST"^^ j 5w{y) QT"^5{x — y), it 
becomes clear that terms not included in Eq. (91) must be 0{gYi,i)^ or in other words 
are small in the weak-coupling limit where SU(N) gauge theory is almost conformal. For 
instance, when neglecting quark masses, bulk viscosity in QCD turns out to be smaller than 
shear viscosity by a factor of gfyM [57]- 

For weakly coupled systems, the form of the non-conformal hydrodynamic equations has 
been investigated in [58] from kinetic theory, but the second-order transport coefficients are 
not known to date. 

For strongly coupled systems, Ref. [59] offers a beautiful example of non-conformal the- 
ories obtained by dimensional reduction of conformal theories. Starting with a conformal 
theory in 2a > D. and reducing to D spacetime dimensions, gives an explicit realization of 
a relativistic hydrodynamic theory where the conformal invariance is (strongly) broken. In 
particular, for this theory the bulk viscosity coefficient C is related to shear shear viscosity 
as 

C .^7/ 1 ,\ 
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and the speed of sound depends on the dimension of the original theory, Cg — The 
relaxation time in the bulk sector tq equals that for the shear sector, 

rn = r^, (109) 
so that one obtains for the limiting velocity Eq. (96) 



max 



+ (110) 



Using the results found in section IV C, . is maximal in the limit of — > 2, where 

?1 , ^ 1. In this limit, vf^"^ — »■ 1, regardless of the value of c?. For D > 2, ^1 , < 1 

and hence vf^'^^ is maximal for the largest possible value of c^, which is = As 
a consequence, one finds that for this class of strongly coupled theories where conformal 
symmetry is broken > 0) the limiting velocity — despite the appearance of Eq. (96) — is 
actually smaller than for a conformal theory in the same number of spacetime dimensions, 
and in particular always smaller than the speed of light. Again, while this does not proof 
that causality is always obeyed in second-order hydrodynamics, it adds to the list of known 
theories where "by coincidence" this turns out to be the case. 

See Ref. [20] for a complete classification of all second-order structures in the energy- 
momentum tensor for non-conformal fiuids. 



V. APPLYING HYDRODYNAMICS TO HIGH ENERGY NUCLEAR 

COLLISIONS 

A. Heavy-Ion Collision Primer 

Relativistic collisions of heavy ions (nuclei with an atomic weight heavier than carbon) 
offer one of the few possibilities to study nuclear matter under extreme conditions in a 
laboratory. The defining parameters for heavy-ion collisions are the center-of-mass collision 
energy per nucleon pair ^/s and the geometry of the colliding nuclei (gold nuclei are typically 
larger than copper, and uranium nuclei are not spherically symmetric). The collisions are 
said to be relativistic once the center-of-mass energy is larger than the rest mass of the 
nuclei, or equivalent ly if y/s/2 is larger than the nucleon mass. For the Lorentz 7 factor of 
the collision, this implies 

^~ md' ~ m ~ 2GeV ' ^ ' 

so typically 7 > 1. Experiments at Brookhaven National Laboratory (AGS, RHIC) and 
CERN (SPS) have provided a wealth of data for Au+Au collisions (AGS, RHIC) and Pb+Pb 
collisions (SPS) ranging in energy from ^/s ~ 2.5 — 4.3 GeV at the AGS over ^/s ~ 8 — 17.3 
GeV at the SPS to ~ 130 - 200 GeV at RHIC. It was found that the number density 
of particles produced in these collisions increases substantially for larger ^/s, indicating a 
similar rise in the energy density [60] , that may allow the study of nuclear matter above the 
deconfincmcnt transition (sec Figure 2). 

For Au+Au collisions at RHIC, the highest energy heavy-ion collisions achieved so far, 
two beams of gold nuclei were accelerated in opposite directions in the RHIC ring and 
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FIG. 2: Left: Particle density (number of particles per unit rapidity, normalized by system size) as 
a function of collision energy (figure from Ref. [60]). Right: The energy density of QCD, calculated 
using lattice gauge theory, shows a strong rise close to the QCD deconfinement temperature (figure 
from [61]). 
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FIG. 3: Schematic view of a relativistic heavy-ion collision. See text for details. 



brought to collide once they reached their design energies. For an energy of a/s = 200 GeV, 
Eq. (Ill) indicates that before the collision the individual gold nuclei are highly Lorentz- 
contracted in the laboratory frame. Thus, rather than picturing the collisions of two spheres, 
one can should think of two "pancake" -like objects approaching and ultimately colliding with 
each other. As a consequence, the duration of the collision itself (which is on the order of 
the nuclear radius divided by the Lorentz gamma factor) is much shorter than the nuclear 
radius divided by the speed of light. Therefore, early after the collision the evolution in the 
directions transverse to the initial beam direction (the "transverse plane") can be assumed 
to be static, and the dynamics is dominated by the longitudinal expansion of the system. 

Being interested in the bulk dynamics of the matter created in a relativistic heavy- 
ion collision, one can divide the evolution into four stages in proper time r = -y/t^ — z^, 
shown schematically in figure 3. Stage I immediately following the collision is the pre- 
equilibrium stage characterized by strong gradients and possibly strong gauge fields [62], 
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where a hydrodynamic description is not applicable. The duration of this stage is unknown 
since the process of equilibration in QCD at realistic values of the coupling is not understood, 
but it is generally assumed to last about 1 fm/c. Stage II is the near-equilibrium regime 
characterized by small gradients where hydrodynamics should be applicable if the local 
temperature is well above the deconfinement transition. This stage lasts about 5 — 10 fm/c, 
until the system becomes too dilute for equilibrium to be maintained and enters stage III, 
the hadron gas regime. The hadron gas is characterized by a comparatively large viscosity 
coefficient [63] , making it ill suited to be described by hydrodynamics, but well approximated 
by kinetic theory [64]. This stage ends when the hadron scattering cross sections become 
too low and particles stop interacting. In stage IV, hadrons then fly on straight lines (free 
streaming) until they reach the detector. 

Assuming the system created by a relativistic collision of two heavy ions becomes nearly 
equilibrated at some instant r = tq in proper time, the subsequent bulk dynamics in stage 
II should be governed by the hydrodynamic equations (23), (91), amended by relevant non- 
conformal terms. To simplify the discussion, in the following these non-conformal terms will 
be neglected, and thus strictly speaking I will not be dealing with real heavy-ion collisions 
but QCD matter in the conformal approximation. However, since the conformal anomaly 
Eq. (107) is small except for a region close to the QCD phase transition [65], there is some 
hope that this approximation does capture most of the important dynamics of real heavy-ion 
collisions. 

To describe the fluid dynamics stage following a heavy-ion collision, one needs to specify 
the value of the hydrodynamic degrees of freedom e^p^u^^n^" at r = Tq, the equation of 
state p — p{e), the transport coefficients 77, r^, Ai^2,3 governing the evolution (91) as well as 
a decoupling procedure to the hadron gas stage at the end of the hydrodynamic evolution. 
None of these are known from first principles, so one typically has to resort to models which 
will be described in the following sections. 



B. Bjorken flow 

The physics of relativistic heavy-ion collisions has been strongly influenced by Bjorken's 
notion of "boost-invariance" [66], or the statement that at a (longitudinal) distance z away 
from the point of (and time t after) the collision, the matter should be moving with a velocity 

— z/t. Neglecting transverse dynamics {v^ — v''^ — 0) and introducing Milne coordinates 
proper time r = \Jt'^ — z"^ and spacetime rapidity ^ = arctanh(2;/i), boost-invariance for 
hydrodynamics simply translates into 

X ^ F .sinh^ zCOsh£ ^ 
^ - , = -u^ -^0 (112) 

T T T 

and as a consequence e^p^u^^n^^ are all independent of ^, and therefore unchanged when 
performing a Lorentz-boost. 

Even though in this highly simplifled model the hydrodynamic degrees of freedom now 
only depend on proper time r, the system dynamics is not entirely trivial. The reason for 
this is that in Milne coordinates, the metric is given by g^y = diag(l, —1,-1, — r^) and hence 
is no longer coordinate- invariant. Indeed, one finds that the Christoffel symbols (IV B) are 
non-zero, 

r|. = ^ , - r (113) 
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so as a consequence the covariant fluid gradients are non- vanishing 

v^u^' = dx + V = r|r«" = ^ ^ ' (11^) 

even though the fluid velocities are constant = (1,0) ! In essence, the Milne coordinate 
system describes a space-time that is expanding one-dimensionally, so that a system at 
rest within these coordinates "feels" gradients from the "stretching" of spacetime, akin to 
the effect of Hubble expansion in cosmology. Unlike in cosmology, however, the spacetime 
described by Milne coordinate is flat, as can be verified by showing that the Ricci scalar 
R — 0. This is important, since one does not want to describe heavy-ion collisions in 
curved spacetime, but rather use the Milne coordinates as a convenient way to implement 
the rapid longitudinal expansion following heavy-ion collisions. Indeed, the covariant fluid 
gradient in Milne coordinates is precisely the same as the one from Bjorken's boost-invariance 
hypothesis, 

V^u^ = - = dt- + dA. (115) 
r r r 

This longitudinal flow (or "Bjorken flow"), together with the assumption = = 0, can 
be seen as a toy model for the hydrodynamic stage following the collision of two inflnitely 
large, homogeneous nuclei. The initial conditions for hydrodynamics at r = tq are then 
completely specified by two numbers: the initial energy density e(To) and one component of 
the viscous stress tensor, e.g. 7r|(ro) (the other components of tt^ are completely determined 
by symmetries as well as u^t:^^ = irjl = 0). For example, in ideal hydrodynamics one finds 
(cf. Eq. (23)) 

De+(e + p)VX = ^re + ^ = (116) 

T 

for the evolution of the energy density (the evolution equations for Du°' are trivially sat- 
isfied). For an equation of state with a constant speed of sound Cg this can be solved 
analytically to give 

e{r) = e{ro) (^y^'^^\ (117) 



r 

Therefore, the energy density is decreasing from its starting value because of the longitudinal 
system expansion, with an exponent that depends on the value of the speed of sound. For 
an ideal gas of relativistic particles = 1/3, giving rise to the behavior e oc t~*^^ that is 
sometimes used in heavy-ion phenomenology. 

Viscous corrections to Eq. (117) have been calculated in first order viscous hydrodynamics 
[67, 68] (the acausality problem discussed in section 11 B does not appear for Bjorken fiow 
due to the trivial fiuid velocities), as well as second order viscous hydrodynamics [27, 69, 72], 
where the equations become 

Moreover, higher order corrections are accessible for known supergravity duals to gauge 
theories [70, 71]. Due to its simplicity, one can expect that Bjorken fiow will continue 
to be used as a toy model of a heavy-ion collisions also in the future, and indeed also I 
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FIG. 4: Rapidity dependence of produced particles (pions), for different collision energies (figure 
from Ref. [60]). Even at highest energies, a plateau at y = (boost- invariance) does not seem to 
emerge. 

will assume rapidity-independence for the remainder of the discussion on the hydrodynamic 
descriptions for simplification. However, it is imperative to recall that experimental data 

by no means supports Bjorken's hypothesis of rapidity independence, as is shown in Fig. 4. 
Rather, the data suggests that the rapidity shape of produced particles is approximately 
Gaussian, independent of the collision energy. This clearly limits the applicability of the 
boost-invariance assumption to the central rapidity region (close to the peak of the Gaussian 
in Fig. 4). 



C. Initial conditions for a hydrodynamic description of heavy- ion collisions 

While pure Bjorken flow assumes matter to be homogeneous and static in the transverse 
(xx = {x,y)) directions, a more realistic model of a heavy- ion collision will have to include 
the dynamics in the transverse plane. This means one has to specify the initial values for the 
hydrodynamic degrees of freedom as a function of x_l. While it is customarily assumed that 
the fluid velocities initially vanish, u^{ro,^±) = 'U^(ro,x_L) = 0, there are two main models 
for the initial energy density profile e(ro,xx): the Glauber and Color- Glass- Condensate 
(CGC) model, respectively. 

The main building block for both models is the charge density of nuclei which can be 
parameterized by the Woods-Saxon potential, 

Pa{x) = — -rr^ (119) 

1 + exp[(|a;| - Ro)/ao\ 

where Rq, gq are the nuclear radius and skin thickness parameter, which for a gold nucleus 
take values of Rq ~ 6.4 fm and Qq ~ 0.54 fm. po is an overall constant that is determined 
by requiring / d^xp^lx) = A, where A is the atomic weight of the nucleus (A ~ 197 for 
gold). In a relativistic nuclear collision, the nuclei appear highly Lorentz contracted in the 
laboratory frame, so it is useful to define the "thickness function" 

/oo 
dzpA{x), (120) 
-oo 
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FIG. 5: Left: Number of participants A'^part in the Glauber model as a function of impact parameter 
b. Right: Spatial eccentricity and area of the overlap region for the Glauber and CGC model, as 
a function of A^part- 



which corresponds to "squeezing" the nucleus charge density into a thin sheet. 

In its simplest version, the Glauber model for the initial energy density profile following 
the collision of two nuclei at an energy ^/s with impact parameter b is then given by 

b b 

e(xx, b) = const x Ta{x + ^^v) ^ Ta{x - -,y) (Tjviv(v^) , (121) 

where (J]^isi{y/s) is the nucleon-nucleon cross section and the constant is freely adjustable 
(see [74] for more complicated versions of the Glauber model). Eq. (121) has the geometric 
interpretation that the energy deposited at position is proportional to the number of 
binary collisions, given by the number of charges at in one nucleus times the number of 
charges at this position in the other nucleus, times the probability that these charges hit 
each other at energy ^/s. Another concept often used in heavy-ion collision literature is the 
number of participants A^part(&) = / '^^2;±^Part(x_L, 6), where 



riPart(x±, b) = nj^art(x±, &) + ^Part(x±, -&) 

/ O-tvatTa (x - 

1- 1 — 



^Part(x±,&) = TA\X + -,y 



A 




;i22) 



Experiments are able to determine the number of participants, but cannot access the impact 
parameter of a heavy-ion collision directly, so the Glauber model A/'part rather than b is 
customarily used to characterize the centrality of a collision (see Figure 5). 

The CGC model is based on the fact that a nucleus consists of quarks and gluons which 
will interact according to the laws of QCD. Accordingly, one expects corrections to the geo- 
metric Glauber model due to the non-linear nature of the QCD interactions. Heuristically, 
one can understand this as follows [75]: At relativistic energies, the nucleus in the labora- 
tory frame is contracted into a sheet, so all the discussion focuses on the dynamics in the 
transverse plane. There, the area vrr^^ of a gluon is related to its momentum Q via the 
uncertainty principle, \rgi\ x \Q\ h = 1, and the cross-section of gluon-gluon scattering in 
QCD is therefore 

a~«,(g2)vrrj^«,(g2)^, (123) 
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FIG. 6: Schematic view of a heavy-ion cohision at impact parameter b in the transverse plane 
(Figure from [73]). 



where is the strong coupUng constant. The total number of gluons can be taken to be 
roughly proportional to the number of partons in a nucleus, and hence also to its atomic 
weight A. Therefore, the density of gluons in the transverse plane is approximately A/ (tcRq), 
where Rq is again the nuclear radius. Gluons will start to interact with each other if the 
scattering probability becomes of order unity. 



A 



A 



a. 



(124) 



which separates 



Therefore, one finds that there is a typical momentum scale Ql 

perturbative phenomena (Q^ ^ Ql) from non-perturbative physics at ^ Ql (some- 
times referred to as "saturation"). The Color- Glass- Condensate was invented [76, 77] to 
include the saturation physics at low momenta Q^ in high energy nuclear collisions. Due 
to the high occupation number at low momenta, this physics turns out to be well approx- 
imated by classical chromodynamics. Despite encouraging progress [78], the problem of 
calculating the energy density distribution in the transverse plane at r = tq using the Color- 
Glass-Condensate has not been solved, the main obstacle being the presence of non-abelian 
plasma instabilities [79, 80]. As a consequence, there only exist phenomenological models 
for the transverse energy distribution in the CGC (which are quite successful in describing 
experimental data, cf. [81]), in particular the model by Ref. [82], which will be referred to 
as CGC model in the following. 

In the CGC model, the transverse energy profile at r = tq is given by 



e(x_|_, h) = const x 



n4/3 



d^XTdY 



[XT,&) 



(125) 



where A^„ is the number of gluons produced in the collision. 
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In order to see the difference between the Glauber and CGC model, one defines the spatial 
eccentricity 



Ah) = '-br-zJr^ (127) 



and overlap area 

S-ovcrlapl^) = 'K^{x^)e{y^)e (128) 

where ()e denote integration over the transverse plane with weight e(x^, h). These quantities 
characterize the shape of the energy density profile in the transverse plane (cf. Fig. 6) and are 
shown in Fig. 5. One finds that the CGC model typically has a larger eccentricity than the 
Glauber model, which will turn out to have consequences for the subsequent hydrodynamic 
evolution. To see this, note that if > 0, the energy density drops more quickly in the 
x-direction than in the y-direction because the overlap region is shaped elliptically. Using 
an equation of state p = p{e) this implies that the mean pressure gradients are unequal, 
dxP > dyp, and according to the hydrodynamic equations (23) one expects a larger fiuid 
velocity to build up in the x-direction than in the y-direction. Since the CGC model has a 
larger Cx than the Glauber model, this anisotropy in the fiuid velocities should be larger for 
the CGC model, as will be verified below. 



D. Numerical solution of hydrodynamic equations 

The hydrodynamic equations are a set of coupled partial differential equations with known 
initial conditions. Typically, it is not known how to find analytic solutions to these set of 
equations, so it is necessary to come up with numerical algorithms capable of solving the 
hydrodynamic equations. As a toy problem, it is useful to study cases where the equations 
simplify, e.g. the assumption of Bjorken fiow discussed in section VB where the hydrody- 
namic equations become a set of ordinary differential equations (118). A standard algorithm 
to solve Eqns. (118) numerically is to discretize time, r = tq+u x At, where tq is the starting 
value, n is an integer, and Ar is the step-size that has to be chosen small enough for the 
algorithm to be accurate, but large enough for the overall computing time to be reasonable. 
With this discretization, derivatives become finite differences, e.g. 

dreir) = , (129) 

and (118) becomes 

n^,r.+i - ^i,n + ^ry + g^^^^^ ^ ^^^^ g^^^ ^ ^^^^ 7r^,„ ^^^^^ ^7r^,„J ^ 

where for simplicity ?7, Ttt, Ai were assumed to be independent of time. Knowing e,p, 7r| at 
step n, the r.h.s. of the above equations are explicitly known (the reason for this was the 
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choice of "forward- differencing" (129)) and hence one can calculate e,p, 7r| at step (n + 1). 
Repetition of this process gives a numerical solution for given stepsize Ar. Since the physical 
solution should be independent from the step size, it is highly recommended to create several 
numerical solutions for different Ar and observe their convergence to a "continuum solution" 
for At 0. 

Unfortunately, the above strategy of simple discretization does not always lead to a well- 
behaved continuum solution. To see this, consider as another toy problem the numerical 
solution f{t,x) to the partial differential equation 



dtf{t,x) = -aod^f{t,x), 



(131) 



where oq is assumed to be constant. Again discretizing time and space as t — to + nAt, x — 
mAx, the derivatives can be approximated by the finite differences 



dtf{t,x)^ 



fn 



At 



dxf{t,x) = 



fn. 



m+1 



fn. 



m—1 



2Ax 



(132) 



which gives rise to the "forward-time, centered-space" or "FTCS" algorithm [83] §19. This 
algorithm is simple, allows explicit integration of the differential equations as in Eq. (131), 
and usually does not work because it is numerically unstable. The instability can be easily 
identified by making a Fourier-mode ansatz for f{t,x) = g*'^"^*-**:'"^^ a^fj calculating the 
dispersion relation lu = uj{k) from the FTCS-discretized Eq. (131) 
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One finds 



Im uj{k) 
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(133) 



(134) 



which signals exponential growth in f{t, x) for all modes k. As a consequence, any numerical 
solution to Eq. (131) using the FTCS algorithm will become unstable after a finite simulation 
time set by the inverse of Eq. (134). 

However, this instability can be cured by choosing a slightly different way of calculating 
the time derivative, namely replacing /„^^ in Eq. (132) by its space average |(/n,m+i + 

/n,m— l) ) 

Q /./, \ fn+l,m fn,m fn,m+l '^fn,m ~l~ fn,m—l fl35l 

This algorithm, known as the "LAX" scheme [83] §19, has a dispersion relation with 



Im u!{k) — 



1 



2At 



In 



COS kAx + [ c^o^— I sm kAx 



(136) 



and hence is numerically stable for Oox^ < 1) e.g. for sufficiently small time steps At. The 



stability of the LAX scheme comes from the presence of the last term in Eq. (135), which 
in "continuum-form" is a second derivative that leads to 



dj{t,x) = -aodJ{t,x) + ^-^dlfit^x) 



(137) 
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FIG. 7: Left to right: snapshots at r = 1,3 and 7 fm/c of the temperature profile T(x,y) for a 
hydrodynamic simulation of a ^/s = 200 GeV Au+Au collision at 6 = 10 fm. The initial spatial 
eccentricity is gradually converted into momentum eccentricity and almost disappears at late times. 



instead of Eq. (131). For sufficiently small Ax, this equation reduces to the original equation, 
so the LAX algorithm indeed converges to the physically interesting solution. But the 
presence of this extra term, which is crucial for the numerical stability, also has a physical 
interpretation: comparing Eq. (137) to the diffusion equation (33) one is led to interpret the 

coefficient as "numerical viscosity". The LAX scheme works where the FTCS scheme 
fails because the viscous term dampens the instabilities, in much the same way that the 
turbulent instability in fluids is damped by the viscous terms [2] §26. Indeed, for ideal fluid 
dynamics numerical viscosity is essential for stabilizing the numerical algorithms. On the 
other hand, viscous fluid dynamics comes with real viscosity inbuilt, so it is tempting to 
conjecture that as long as or ^ are finite and At is sufficiently small, numerical viscosity 
is not needed to stabilize the numerical algorithm for solving the hydrodynamic equations, 
and the simple FTCS scheme can be used. Indeed, at least for the problem of heavy-ion 
collision, this strategy leads to a stable algorithm [84-87]. 

Aiming to solve the hydrodynamic equations in the transverse plane (assuming boost- 
invariance in the longitudinal direction), one first has to choose a set of independent hy- 
drodynamic degrees of freedom, e.g., e,u^,uy,7i^^,n^y,n^y for which initial conditions are 
provided along the lines of section VC. Only time derivatives to first order of these six 
quantities enter the coupled partial differential equations (23), (91), so that formally one can 
write the hydrodynamic equations in matrix form 



A 



drU^ 



(138) 



where A, b are a matrix and vector with coefficients that do not involve time derivatives. 
Using the FTCS scheme to discretize derivatives, and matrix inversion to solve (138), the 
value of the independent hydrodynamic degrees of freedom at the next time step are explicitly 
given in terms of known quantities (once the equation of state and hydrodynamic transport 
coefficients are specified). Reconstructing all hydrodynamic fields from the independent 
components and repetition of the above procedure then leads to a numerical solution for 
the hydrodynamic evolution of a heavy-ion collision for given Ar, Ax as long as rj > (in 
practice, values as low as rj/s ~ 10~^ are stable with reasonable Ar). The convergence 
of these numerical solutions to the continuum limit is explicitly observed when choosing 
a series of sufficiently small step sizes Ar, Ax. Snapshots of the temperature profile in a 
typical simulation are shown in Fig. 7. 
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FIG. 8: Time evolution of the spatial and momentum anisotropies, Eq. (127) and Eq. (139), 
respectively, in the Glauber and CGC model for a ^/s = 200 GeV Au+Au collision at 6 = 7 fm 
(Figure from [72]). 



Fig. 7 also displays the gradual reduction of the eccentricity (the shape of the temperature 
profile in the transverse plane becomes more and more circular as time progresses). The 
eccentricity corresponds to a spatial anisotropy in the pressure gradients which is converted 
by hydrodynamics into a momentum anisotropy (fluid velocities ^ u^). In analogy to the 
definition of the spatial eccentricity (127), it is therefore useful to introduce the concept of 
momentum anisotropy 

The time evolution of the eccentricity and momentum anisotropy in the Glauber and CGC 
model are shown in Fig. 8. As discussed in section VC, the higher initial eccentricity in the 
CGC model gets converted in a larger momentum anisotropy. 



E. Preeze-out 



Experiments in relativistic nuclear collisions measure momentum distributions of particles 
(pions, kaons, protons, etc.), whereas hydrodynamics deals with energy density, pressure 
and fluid velocities. Clearly, in order to make contact with experiment, the hydrodynamic 
degrees of freedom need to be converted into measurable quantities, which is often called 
the "freeze-out". The connection between hydrodynamics and particle degrees of freedom 
is provided by kinetic theory, which was discussed in section III. In particular, one requires 
the hydrodynamic and kinetic theory energy momentum tensor at freeze-out to be the same, 

^kinetic theory = / dxP^P^fif, t, x) = T^^^^^ , (140) 

where for small departures from equilibrium the explicit form of / in terms of hydrodynamic 
degrees of freedom is provided by Eq. (59). Once /(p, t,a;) is known, one can construct the 
particle current from kinetic theory 



= / dxp^f{p,t,x) , 



(141) 
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which will be used to construct particle spectra that can be compared to experimental 
measurements. 

Freeze-out from hydrodynamic to particle degrees of freedom is expected to occur when 
the interactions are no longer strong enough to keep the system close to thermal equilibrium. 
Below the QCD phase transition, this happens, e.g., when the system cools and viscosity 
increases [63] until the viscous corrections in (59) become too large and a fluid dynamic 
description is no longer valid. In practice, this is hard to implement, so simplifled approaches 
such as isochronous and isothermal freeze-out are often used (see, however, [88, 89]). All of 
these have in common that they define a three dimensional hypersurface E with a normal 
vector dT,'^ such that the total number of particles after freeze-out is given by the particle 
current (141) leaving the hypersurface E'*, 

N^J n^dE^ = J dxfip, t, x)p^dJ:^ . (142) 

For energy densities sufficiently below the QCD phase transition, the energy momentum 
tensor is well approximated by a non-interacting hadron resonance gas [92]. This translates 
to 

dx = E(2^^ + 1)(2^^ + ^)-0j3^(P'P^ - ^i)np'), (143) 

where the sum is over all known hadron resonances [93] and Sj, gi are the spin and isospin 
of a resonance with mass m^. As a consequence. 



N^^I^^^Ui^] (144) 



where 

odN\ 



p 



I d^^pi^fip, t, x), di^ {2s, + l){2gi + 1) , (145) 



is the single-particle spectrum for the resonance i. Eq. (145) is the generalization of the 
"Cooper- Frye freeze-out prescription" [94] to viscous fiuids with / given by Eq. (65). 

Arguably the simplest model is isochronous freeze-out, where the system is assumed to 
convert to particles at a given constant time (or proper time). While fairly unrealistic, 
it allows a rather intuitive introduction of the general freeze-out formalism: constant time 
defines S'^(t, x, y, z) in the hydrodynamic evolution which is parametrized by t = const. The 
normal vector dH'^ on this hypersurface is given by [90, 91] 

dT.'^d^^dTP , , , 
= e^apy-Q^-Q;^-Q^dxdydz , (146) 

where Sf^cep-y is the totally antisymmetric tensor in four dimensions with 60123 = +1- A simple 
calculation gives rfS^ = {l,0)d'^x and therefore the momentum particle spectra are easily 
obtained by integration of the distribution function, 

p'^dE^ = p'd'x , (^^^ = (^ / ^ = ^) • 
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Slightly more realistic is isochronous freeze-out in proper time, r = const, where the 
freeze-out surface = (r cosh ^, x, y, r sinh ^) is parametrized by x, y, ^, because this incor- 
porates Bjorken flow. Introducing rapidity Y = arctanh(p^/p°) and m± = ^ m? -\- p\ for 
convenience, a short calculation for the normal vector dS^ gives [72] 

dH^p^" = Tmi_ cosh (y - OdxdydC . (148) 

Considering for illustration a Boltzmann gas in equilibrium with /eq = exp [— p'^u^/T], for 
vanishing spatial fluid velocities one has 

(^) ^ (^p™"^ / ^^^y^^ ~ ^) ^-^P [-'m± cosh (y - ^)/T] 

= (149) 

while for non- vanishing fluid velocities with azimuthal symmetry {u^{r) = u^ir) — 
u^{r)/V2) a short calculation gives [84] 



d^pj^ (27r) 



r-./^-X.(^)/o ^ , (150) 



where K{z),I{z) are modified Bessel functions and the transverse radius r = y/x'^ + y'^ 
has been introduced for convenience. Comparing the integrands in Eq. (149), (150) when 
transverse momenta p± — \p_l\ are much larger than the temperature T or mass m, one 
finds 

Ki{m^/T) exp [{m±{l - u^) + p±u')/T] ^ 

Ki{m±u^ /T)Io{p±W /T) ^ u^u'-p^/T ^ ' 

~ 0{1). This means that the presence of u'' > 0, or "radial fiow", leads to particle 
spectra which are "fiatter" at large pt- This is confirmed by numerical simulations [95]. 

For a Boltzmann gas out of equilibrium and Bjorken fiow only {u^ — l,u^ — the viscous 
correction to the distribution function (65) is 

tt'^^p^ _ (tt^^ + T^yy)pi + 27r««mi/T2 sinh^ {Y - _ T^f (pi - 2m\ sinh^ (y - S,)) 



2(e + p)T2 4(e + p)T2 A{e + p)T^ 

so that the single particle spectrum becomes [27] 



JdN\ 2d, f. 



\Tj « 4(e + p)r2 



(152) 



(153) 



Since for Bjorken flow typically 7r| > 0, this implies that viscous corrections tend to have 
the same effect of making particle spectra flatter at large Pt, which hints at the difficulty of 
extracting viscosity and radial flow from experimental data [85] . More information is needed 
to disentangle these effects, so one decomposes the particle spectra into a Fourier series with 
respect to the azimuthal angle in momentum (f)p — arctan (p^/p^) [96], 

(^°^) ^ M\P±\,b) [I + 2v2{\p±\,b) cos2(l)p + 2vi{\p^\,b) cosicPp + . . .] , (154) 
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FIG. 9: Entropy per particle as a function of temperature for a gas of a realistic set of massive 
resonances with Bose-Einstein/Fermi-Dirac statistics (dotted line). Also shown are results when 
approximating by Boltzmann statistics (dashed line) and total initial entropy over final multiplicity 
from a numerical hydrodynamics simulation with isothermal freeze-out using Boltzmann statistics 
[72] (squares). The numerical results prior to the decay of unstable resonances is in fair agreement 
with the analytic prediction. The decay of unstable resonances produces additional particles, 
leading to a smaller ratio of initial entropy to final multiplicity (circles). 



where the coefficients v^, f4 are referred to as "elliptic" and "hexadecupole" flow [97], respec- 
tively. V2iVi^ and even higher harmonics were measured experimentally for a/s = 200 GeV 
Au+Au collisions [98] and may be useful to distinguish between flow and viscous effects. 

A more realistic criterion than isochronous freeze-out is to assume decoupling at a 
predefined temperature (isothermal freeze-out). In this case the hypersurface S can be 
parametrized by ^, = arctan| and a time-like coordinate < c < 1 with c = cor- 
responding to the center x = y = of the transverse plane. Assuming boost-invariance 
in the longitudinal direction, this leads to Y)^ = S^(r(c) cosh ^, x(c, 0), ?/(c, 0), r(c) sinh^), 
and the normal vector is evaluated analogous to Eq. (146) (cf.[72]). The resulting single 
particle spectra are then given by Eq. (145), where it may be convenient to change variables 
c = ^ in the integral 



dc=- , (155) 

Jto Tfo — To 

where tq, Tfo correspond to the start and end of the hydrodynamic evolution. For isothermal 
freeze-out at a temperature Tfo, kinetic theory specifies the entropy density s = ^ and 
the number density n = n^u^ of particles. In particular, for a massive Boltzmann gas in 
equilibrium Eq. (57), (141) lead to 

^ (23, + l)(2(7, + l) 3 /mA ^ (2g, + l)(2(7, + l) ^ (m,\ 

which can be used to quantify s/n, the amount of entropy each resonance degree of freedom 
is carrying. For extremely high temperatures s/n 4, which is the known limit for a gas 
of massless relativistic particles [99], but for temperatures below the QCD phase transition 
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and a realistic set of hadron resonances [93], s/n depends on Tfo (see Fig. 9). For ideal 
hydrodynamics the total entropy S in the fluid is conserved (27), and hence the total number 
of particles created by an isothermal freeze-out should be given by iV = ^5, which provides 
a non-trivial check on numerical codes. 

After freeze-out, the hadron gas dynamics may be described by a hadron cascade code 
such as [64]. A more simplistic approach is to assume that particles stop interacting after 
freeze-out, but unstable particles are allowed to decay, which changes the spectra of stable 
particles [100, 101]. The decay of unstable resonances can be simulated using pubhc codes 
such as [102] and leads to particle production, as can be seen in Fig. 9. 

F. Viscous effects and open problems 

Ideal hydrodynamic simulations have been used quite successfully in the past to describe 
the properties of the particle spectra produced in relativistic heavy-ion collisions [60, 103- 
105] (for reviews, see e.g. [95, 106]). Viscous effects have only been studied more recently: 
the presence of viscosity leads to viscous entropy production given by Eq. (27), which in- 
creases the total multiplicity for fixed initial entropy. The amount of viscous entropy pro- 
duction depends on the hydrodynamic initialization time Tq [82], and for Tq ~ 1 fm/c is on 
the order of 10 percent for rj/s = 0.08 [85, 108]. 

Viscosity also leads to stronger radial flow, which increases the mean transverse momen- 
tum of particles [84, 85, 109]. Maybe more importantly, the presence of shear viscosity 
strongly decreases the elliptic flow coefficient V2. After some initial disagreement, several 
different groups now agree on the quantitative suppression of elliptic flow by shear viscosity, 
as is demonstrated in Fig. 10. Unfortunately, this does not directly constrain the 77/5 of 
hot QCD matter because the overall size of elliptic flow (which is proportional to the final 
momentum anisotropy Cp) is dictated by the initial spatial eccentricity, which is unknown 
(cf. Fig. 8, see also Ref. [72]). 

Many open problems remain, such as 

• Exploring the effects of bulk viscosity. One study [110] suggests that bulk viscosity 
may become large close to the QCD phase transition (however, see [51, 56, 111]), 
which would have important consequences for the hydrodynamic evolution. First 
phenomenological steps in this direction have been taken in [112], but it would be 
worthwile to have a classification of all the non-conformal terms (including an estimate 
of their importance/size) in the hydrodynamic equations. 

• Implementing finite baryon chemical potential in the viscous hydrodynamic simula- 
tions. On the one hand needed to describe the asymmetry in the baryon/anti-baryon 
multiplicities, the viscous hydrodynamic evolution in the vicinity of a possible QCD 
critical point could on the other hand help to guide experimental searches for this 
critical point (cf.[113]). 

• Combining viscous hydrodynamics with a hadron cascade code to more realistically 
describe the freeze-out process. Ideally such a hybrid code would make the choice of a 
freeze-out temperature superfluous, eliminating one model parameter (cf. [114-1 16]). 

These (and other) problems are straightforward to solve in the sense that they do not 
require fundamentally new ideas, but "only" hard work. However, there are also at least 
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FIG. 10: Reduction of elliptic flow coefficient due to shear viscosity: different groups agree that at 
■p^ = 1.5 GeV, there is a ~ 20 percent reduction of W2 for rj/s = 0.08 (Figures from [86, 89, 107, 108], 
clockwise from upper left.) 



two problems outside the framework of hydrodynamics, which would have to be solved in 
order to claim a complete understanding of high energy nuclear collisions: 

• What is the value of the initial spatial eccentricity in high energy nuclear collision? 
Can it be calculated or measured without having to rely on models (Glauber/CGC)? 
Since the spatial eccentricity controls the amount of elliptic flow generated in a rela- 
tivistic nuclear collision, knowing Cx seems necessary to quantify the viscosity of hot 
nuclear matter. 

• How and when does the system equilibrate? An answer to this question would give 
a well defined value to the hydrodynamic starting time tq as well as the eccentricity 
at this time. Currently both tq and exiro) are "guessed", with no solid arguments for 
any particular value. 

Nevertheless, the striking ability of viscous hydrodynamics to describe the momentum 
spectra of the majority of particles, including the elliptic flow coefficient, in the highest 
energy Au+Au collisions at RHIC (see Fig. 11) make relativistic nuclear collisions an ideal 
application for the old, new, and future developments in hydrodynamic theory. 
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FIG. 11: Centrality dependence of particle spectra (multiplicity, mean momentum and elliptic 
flow) as well as momentum dependence of charge hadron elliptic flow, for a viscous hydrodynamic 
simulation of a -^s = 200 GeV Au+Au collision using the CGC model compared to experimental 
data [119-121] (Figures from [72]). 



VI. CONCLUSIONS 

Relativistic viscous hydrodynamics is an effective theory for the long-wavelength behavior 
of matter. The relativistic Navier-Stokes equation would do justice to this long-wavelength 
behavior, but does not lend itself easily to direct numerical simulations because of its ul- 
traviolet behavior. Generalizations of the Navier-Stokes equation including second-order 
gradients have been proposed to cure this difficulty, and indeed provide a phenomenological 
regularization of the Navier-Stokes equation that can be simulated numerically unless the 
regularization parameter is too small. 

Interestingly, performing a complete gradient expansion to second order reproduces this 
attractive feature of regularizing the Navier-Stokes equation, besides having the benefit 
of constituting an improved approximation of the underlying quantum field theory. For all 
theories where the regularization parameter obtained from this gradient expansion is known, 
its value is such that the ultraviolet behavior of the second-order hydrodynamic equations 
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is benign. It is not known whether this is a coincidence. 

The second order hydrodynamic equations have been apphed to the problem of high 
energy nuclear collisions, offering a good description of the experimentally measured particle 
spectra at low momenta. Further work is needed to extract material parameters of hot 
nuclear matter, such as the shear viscosity coefficient, from experimental data. 

Many other applications of second order hydrodynamics are possible, e.g. in astrophysics 
(viscous damping of neutron star r-modes [117]) or cosmology (effects of bulk viscosity [118]). 

Whether in the formulation I have described in these pages, or not, one thing is certain: 
relativistic viscous hydrodynamics is here to stay. 
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APPENDIX A: PROOF OF CAUSALITY OF MAXWELL-CATTANEO TYPE 

EQUATIONS 



Let us first establish that the diffusion-type equation 



df5uy 



eo+Po 



d'juy^f(t,x), 



(Al) 



which was discussed in section II B for the homogeneous case f{t,x) = 0, violates causality. 
To this end, let us calculate the retarded Green's function G(x, x'), x = (t, x), of the 
differential operator (Al), 



G(x,x') = 5^(x-x'), 



Doing a Fourier-transform of G one finds 



G(x,x') 



^2j^ ^—iu>(t—t')+ik{x—x') 



(27r) 



-ioj + vk^ 



(A2) 



(A3) 



which can be solved by the usual contour integration methods and Gaussian integration, 

{x — x')^ 



G(x, x') = I ^ ■ = exp 



The solution to Eq. (Al) is then 



5v}>{t,x) = / ciVG'(x,x')/(x') = -^i2=exp 



X 



Avt 



(A4) 
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where the system was started with an initial "kick", f{t,x) = S{t)S{x). One can see that 
for any finite time t > 0, the perturbation is non-vanishing for all values of x, not only for 
those X < t. This obviously violates causality. 

Considering instead of the diffusion- type equation (Al) the Maxwell-Cattaneo law 



dtSu'' + — — 

eo+Po 

the Green's function has to fulfill 
and hence is given by 



G'(x,x') = -5\:>c-k') 



^(x, x') 



-iLo{t—t')+ik{x—x') 



The frequency integration proceeds as before, and one finds 



~+»^ dk i 

-oo+ie 27r Ttt 



^ik{x—x') 



-iuj+{t-t') 



-iw-(t-t') 



UJ~ 



(A5) 



(A6) 



(A7) 



(A8) 



where 2r7ra;^ = — i ± -\/4T7ri^A;^ — 1. The integral over A; is chosen in the upper half-plane and 
the branch cut of the square- root is chosen to run from k — — (4t^z/)~^/^ to A; = (4r^i/)~^/^ 
along the real axis [122] §7.4. To evaluate the integral 



/. = e 



oo+ie I exp 



— y) 



: (t-t') 



-oo+ie 27r 



note that \ix > x' + ^^:^a/4tv^, then the contour can be closed by a semicircle in the upper 
half plane, giving a vanishing contribution since there are no singularities in that half plane. 
For X < x' + ^^\/^ Tt^v on the other hand, the contribution will not vanish. It can be 
calculated by using a table of Laplace transforms [122]§7.4, [123]§11, giving 



— (X — X 



{t - t'f {x - x') 



i\2 



(A9) 



where Iq{x) is a modified Bessel function. Similarly, one can calculate the other component 
of Eq. (A8), so that one finds 



G'(x,x') = d{t-f)d 



\t-t'fv 



{x — x) 



'\2 



t-t' 



{t - t'f {x - x'f 



V4l^ y\ 4r2 4i/T^ 



. (AlO) 



Prom the step-function in Eq. (AlO), one can easily convince oneself that the solution 
5uy{t,x) to Eq. (A5) is confined to \x\ < tv^^, where — ^ ^ (eo+po) " coincides with 
the limit found in Eq. (38). The difference between the Maxwell-Cattaneo solution and the 
diffusion equation is highlighted in Fig. 12, where \/47n4G(x, 0) is plotted for t = lOz/ as a 
function of x. One can see that G(x, 0) has non- vanishing support in the region excluded by 
causality for the diffusion equation, while this does not happen for the Maxwell-Cattaneo 
law. 
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FIG. 12: Green's function for the diffusion equation and Maxwell-Cattaneo law for r^r = i'. See 
text for details 

APPENDIX B: NOTATIONS AND CONVENTIONS 

This appendix is a collection of notations and conventions used in the main part of the 
article. 

• The metric sign convention is (+,—,—,—) 

• Projectors: 

A^^ = g^'^ - u^u'^ , Pi^l = K^^p + ^^pK - l^^^^^p , (Bl) 

with properties Uf.A'"' = u^P^^p = 0, g^uPap = 0. 

• Derivatives: 

D = u^'D^, V^ = A^D^, D^ = u^D + V^, (B2) 

where is the geometric covariant derivative that reduces to for flat space. 

In the non-relativistic, flat-space limit, 

D = dt + v-d + 0{\v\^) , W = -d + 0{\v\), (B3) 

which supports the interpretation of time-, and space-like derivatives for D and V, 
respectively. 

• Brackets: 

A«^B^> = Pl^^l^At'B'' , (B4) 

which are used to define e.g. the vorticity, VL^p = ViaUp]. Note that the above 
definition of A'^'^B^^ differs from others (e.g. [29]) by a factor of 2. 
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